Metabolic Network Modelling: Including Stochastic Effects 
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We propose to model the dynamics of metabolic networks from a systems biology point of view 
by four dynamical structure elements: potential function, transverse matrix, degradation matrix, 
and stochastic force. These four elements are balanced to determine the network dynamics, which 
gives arise to a special stochastic differential equation supplemented by a relationship between the 
stochastic force and the degradation matrix. Important network behaviors can be obtained from 
the potential function without explicitly solving for the time-dependent solution. The existence of 
such a potential function suggests a global optimization principle, and the existence stochastic force 
corresponds natural to the hierarchical structure in metabolic networks. We provide theoretical evi- 
dences to justify our proposal by discussing its connections to others large-scale biochemical systems 
approaches, such as the network thermodynamics theory, biochemical systems theory, metabolic con- 
trol analysis, and flux balance analysis. Experimental data displaying stochasticity are also pointed 
out. 



I. INTRODUCTION 

The successful completion of Human Genome Project 
reveals that to understand the molecular base of diseases 
a huge amount of information is needed from every level 
of description: from DNA, protein, to function^. It be- 
comes evident that expertise beyond traditional biology 
is needed to decipher the genomic message. In metabolic 
and physiological studies, there has already been an ac- 
tive tradition of interaction among biologists, chemists, 
physicists, engineers, and computer scientists^. Those 
large scale integrations call for a rethinking on the prac- 
tice of biology. The emergence of systems biology is the 
response to such demand^. 

The systems biology study of metabolic networks may 
be decomposed to be four mutually interconnected parts: 

1) . Biological experiments. It may be in response to a de- 
mand for finding a cure for a particular disease or finding 
a method for health care at the preventive level. Both la- 
bor intensive and optimized designed experiments would 
be carried out. It may be an industrial demand to in- 
crease the yield of a product, such as to increase the pro- 
duction of desired biomass from methanol^. The better 
metabolic network and its related gene regulatory net- 
work needs to be understood before a proper genetic and 
metabolic engineering. Or, the goal is more on biology: 
we simply wish to understand the biological principles in 
terms of molecular, physical, and chemical mechanisms 
behind a given biological process^. To address those 
questions, new tools, both experimental and computa- 
tional, as well as new theoretical framework, are needed. 

2) . Novel experimental technologies. Thus, there is an 
urgent demand for better technologies, particularly in the 
directions of high throughput, real time, and in vivo at 
cellular level. Those technologies should provide us new 
interrogative power to understand biological mechanism, 
with new and more data. 

3) . Biocomputation and bioinformatics. The new and 
large scale data sets call for their proper annotations and 
presentations, the task of biocomputation and bioinfor- 



matics. There have been tremendous efforts recently in 
this direction. 

4). Conceptual and mathematical frameworks. Hence, 
all those new developments call for more powerful concep- 
tual and mathematical frameworks to facilitate a better 
communications among experts of theoretical and exper- 
imental biologists and to help achieve the goals. 
The purpose of the present article is to address the ques- 
tion of mathematical framework from a biophysics and 
biochemistry point of view, motivated by our study in 
systems biology, particularly in the theoretical study of a 
gene regulatory network^ and in the experimental study 
of a metabolic network 4 . 

In addition to large scale biological data sets, more 
time sequence data have been obtaining and in vivo data 
plays increasingly important roles. They often show that 
data are stochastic in natur o 7 i 8 i 9 , not due to our inabil- 
ity to obtain accurate data. They carry molecular and 
cellular biological information. We further notice that 
the biological processes are physical and chemical. This 
implies that we must work with constraints, such mass 
conservation on biomass, the energy conservation, as well 
as other thermodynamic constraints. Hence there exist 
several requirements for the modelling: We demand it to 
provide insights and understandings on bio-structure, to 
perform exploring simulations, to interpret and evaluate 
of measured data, to make prediction and help design 
further experiments, and to optimize the whole process 
to achieve the goal. For those ambitious objectives, we 
believe a proper mathematical framework may be based 
on the stochastic differential equations. This is a mid- 
dle level description: It can be linked to both the master 
equation type description with explicit discrete nature of 
chemical processes and the Fokker-Planck equation type 
with continuous variables in both time and state space. 
One advantage of present description is that a very spe- 
cific form can be spelt out: the four dynamical elements 
formulation of stochastic differential equations. They can 
be directly connected to experimental data. The role 
played by stochascitity is emphasized in this description. 
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In section II a discussion of the dynamical components 
will be given. In section III the specific form will be given 
and its connection to more classical form of stochastic 
differential equation in chemistry is obtained. In section 
IV a further discussion of dynamical elements in sight 
of quantitative equations will be given. Applications 
are discussed in section V. In section VI the connection 
of present framework to several successful approaches 
in metabolic network modelling, the network thermo- 
dynamics theory, biochemical systems theory, metabolic 
control analysis, and flux balance analysis. We conclude 
in section VII. 



II. FOUR DYNAMICAL COMPONENTS IN 
METABOLIC NETWORKS 

Because biological processes are based on physics and 
chemistry, we would naturally assume that they are con- 
tinuous in time. Further, because we are focused on the 
metabolic networks, we would like to assume that the 
most important ingredient in the mathematical descrip- 
tion would be a set of coupled chemical reaction rate 
equations, which may be presented by a set of stochastic 
differential equations. In fact, this has been the working 
assumption underlying all the mathematical metabolic 
network modelling 1 —. This middle level description can 
link to other descriptions of either the complete dis- 
crete master equation type or of the complete continuous 
Fokker-Planck equation type. The important question is 
whether or not a more useful mathematical structure ex- 
ists. 

To address the dynamical structure question of a 
given network describing the macromolecular synthesis 
in metabolic pathways, we may intuitively classify its dy- 
namical components into four basic elements according 
to their roles in the providing or spending of the resource 
in a metabolic network: the driving force, the transverse 
force, the degradation, and the noise. There has already 
been lots of previous attempts in this direction in the con- 
tent of non-equilibrium thermodynamics^. The driving 
force is the active part in such macromolecular synthesis, 
which would include the external force and the internal 
structure determined by the metabolic network architec- 
ture. It provides the necessary resource, such as the en- 
ergy, the mass flow, or other needed chemical ingredients, 
in order to maintain the desired function of the network. 
It is closely related to the robustness of the network, act- 
ing as a backbone or skeleton of network dynamics. It 
is the gradient of the potential function of the network. 
The transverse force describes the process of relocation 
of the resource from one part of the network to another, 
and is responsible for phenomena such as the delay re- 
sponse and the oscillation. The effects of both driving 
and transverse forces alone are reversible: no change of 
the potential function of the network. The dissipation 
describes how the resource is consumed during the net- 
work dynamical process, including the consumption of 



metabolites, dissipation of energy, and degradation of 
macromolecules. It is an irreversible process, and can- 
not be avoided for a complex metabolic network. Owing 
to dissipation, the network cannot maintain its designed 
function if external and internal resources are finished. 
Finally, for any complex network, there always exist fac- 
tors that cannot be known for certain or whose detailed 
dynamical behavior is not our primary interest. Those 
factors form the fourth element: the noise or the stochas- 
tic force. Both dissipation and noise give the network the 
ability to adapt to the optimal state. We will return for 
further discussion of those four elements in Section IV. 

Interestingly, a complex network has only these four 
dynamical elements within present mathematical frame- 
work. Their interplaying gives rise to rich behaviors 
which we will explore elsewhere. In the present article, 
we will present the underlying quantitative mathemati- 
cal framework for the dynamical structure in terms of two 
generic dynamical principles: The equivalence principle 
that all the four elements must be balanced to describe 
the network dynamics; The stochasticity-dissipation rela- 
tion that a constraint exists on the stochastic force. Such 
a constraint implies that the stochastic force must be the 
integral part of the network dynamics, a view most evi- 
dent in biology^. The explicit identification of the four 
dynamical elements of a network leads to a better under- 
standing of the network behavior: The final stationary 
distribution can be directly read out from the driving 
force; the time scales in the network becomes explicit; 
and a long sought optimization principle for a network is 
immediately suggested. 



III. STOCHASTIC DIFFERENTIAL 
EQUATIONS TO DESCRIBE METABOLIC 
NETWORK DYNAMICS 

The questions naturally arise that how do we describe 
the network dynamics quantitatively and how do we iden- 
tify the four elements? To be specific, let us consider 
an n component network. The n components may be 
the relevant the numbers of metabolites in a metabolic 
network or the relevant proteins in a gene regulatory net- 
work pathway- or other quantities specifying the network 
The value of j th component is denoted by qj. The n 
dimensional vector q r = (qi, q^ q n ) is the state vari- 
able of the network. Here the superscript r denotes the 
transpose. Let fj (q) be the deterministic nonlinear force 
on the j th component, which includes both the effects 
from other components and itself, and Q (q, t) the ran- 
dom force. For simplicity we will assume that fj is a 
smooth function explicitly independent of time. The net- 
work dynamics may be generally modelled by the set of 
stochastic differential equations, the stochastic chemical 
rate equation a 10 ! 11 : 

ft = /i(q) + 0(<i*). (i) 
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with j — 1,2, ...,n. Here <jj = dqj/dt. Without loss of 
generality the noise will be assumed to be Gaussian and 
white with the variance, 

(Ci(q,t)0(q,0) = 2A i (q)^-t') ! (2) 

and zero mean, ((j) = 0. Here S(t) is the Dirac delta 
function, and (...) indicates the average with respect to 
the dynamics of the stochastic force. By the physics 
and chemistry convention the semi-positive definite sym- 
metric matrix D = {Dij} with i,j = l,2,...,n is the 
diffusion matrix. Eq.(l) is the standard stochastic dif- 
ferential equation, or the standard stochastic differential 
equation in physics and chemistry. If an average over the 
stochastic force (, a Wiener noise, is performed, Eq.(l) is 
reduced to the deterministic equation in dynamical sys- 
terns, (q) = (f(q)) = f«q)). 

With Eq.(l) and (2) as our starting point to identify 
the dynamical elements, three immediate problems need 
to be addressed. First, the stochastic force £ appears to 
be simply introduced by hand. It is not obvious how it 
can be related to the deterministic force f . Second and 
more serious, apart from the apparent ad hoc stochastic 
force, in Eq.(l) there are no sign of other three dynamical 
elements: the driving force, the transverse force, and the 
dissipation. Third and worse, both the divergence and 
the skew matrix of the force f are in general finite: 

9-f^O, 9xf^0. (3) 

Here the divergence is explicitly d ■ f = Y^j=i dfj/^Qj = 
tr(F), and the skew matrix is twice the anti-symmetric 
part of the force matrix F: (d x f)y = Fji — Fij with 
Fij = dfi/dqj , i,j = l,2,...,n. The finiteness of the 
divergence leads to that the phase space volume is not 
conserved. Dissipation is hence implied. The finiteness 
of the skew matrix has been characterized as the hall- 
mark of the network in a state far away from thermal 
equilibrium^. 

We propose, and have demonstrated elsewhere^, that 
there exists a decomposition such that Eq.(l) can be 
transformed into the following form: 

[5(q)+T(q)]q=-au(q} + e(q > t) I (4) 

with the semi-positive definite symmetric matrix S, the 
anti-symmetric matrix T, the single-valued scalar func- 
tion u, and the stochastic force £. Here d is the gradient 
operator in the state variable space. It is straightforward 
to verify that the symmetric matrix term is 'degradation': 

q r S(q)q > ; 

and the anti-symmetric part does no 'work': 

q T T(q)q = , 

therefore non-decaying. Hence, it is natural to identify 
that the degradation is represented by the symmetric ma- 
trix S, the friction matrix, and the transverse force by 



the anti-symmetric matrix T, the transverse matrix. The 
driving force is clearly represented by the gradient of the 
scalar function u, which acquires now the meaning of po- 
tential energy and will be named the network function. 
There are indeed four dynamical elements for the network 
dynamics. Eq.(4) states that they must be balanced in 
order to produce the network dynamics. 

Eq.(4) is our key for the identification of dynamical 
structure. It is the form of stochastic differentiation 
equation which we believe would be suitable to directly 
analyze metabolic networks. Nevertheless, without fur- 
ther constraint, the four dynamical elements, if exists, are 
not unique. This may be illustrated by a simple count- 
ing: There are four apparent independent quantities in 
Eq.(4), while there are only two in Eq.(l). For this rea- 
son the decomposition from Eq.(l) to (4) will be called 
the singular decomposition. In order to have a unique 
identification, we choose to impose the constraint on the 
stochastic force and the symmetric matrix in Eq.(4) in 
the following manner: 

<e(q,*)r(q,0>=2S(q)^-i'), (5) 

and (£(q)) = 0. This constraint is consistent with the 
Gaussian and white noise assumption for £ in Eq.(l), 
and will be called the stochasticity-dissipation relation. 
The constrained singular decomposition will be called the 
gauged singular decomposition. Eq. (4) and (5) will be 
called the normal stochastic differential equation, or the 
normal stochastic differential equation. 

Two comments are in order. First, it is evident that 
any matrix can be decomposed into a symmetric and an 
antisymmetric matrices. What is unique in Eq.(4) is that 
the symmetric matrix in Eq.(4) has to be semi-positive 
definite, as guaranteed by Eq.(5). Therefore, the matrix 
[S + T] at left hand side of Eq.(4) is not an arbitrary 
matrix. It has a built-in structure which makes the net- 
work tend to the minimum of potential function u(q), 
that is, the tendency for optimization. Since the poten- 
tial function u(q) plays the same role as energy function 
in physics and chemistry, the typical optimization pro- 
cedure, such as the simulation annealing, may be used 
here to find the global minimum of the network potential 
function. This leads to the second comment that with- 
out the stochastic effect in Eq.(l), no unique potential 
function in Eq.(4) can be determined. This implies that 
no global optimization can be found in the absence of 
stochastic effect. 

The connection from Eq.(l), the typical form of chem- 
ical rate equation, to Eq.(4), the form of the present pro- 
posal, may be expressed in following equations^: 

u = - f dq'.[G-Hq')f(q')] 
S = [G- 1 (q) + (G T )- 1 (q)]/2 . (6) 
T = [G-\^-{Gr)-^]/2 

Here the 

Here is the auxiliary matrix function G(q) = [5(q) + 
T(q)] -1 is the solution of following equations: 

d x [G- X f(q)] = , (7) 
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and 

G + G T = 2D . (8) 

With the network function it(q) similar to a potential 
energy, the stationary distribution function po(q) for the 
state variable is expected to be a type of Boltzmann- 
Gibbs distribution: 

Mq) = \ cxp{-u(q)} , (9) 

with the partition function Z — J d n q exp{— u(q)}. This 
is one of most useful results from Eq. (4) : It allows a direct 
comparison to stochastic experimental data at steady 
states. 



IV. FURTHER CLARIFICATION 

With above clear identification of four dynamical ele- 
ments in a complex network from Eq. (4) and the demon- 
stration of its equivalence to the tradition formulation 
of Eq.(l), a further discussion of them is given in this 
section. 

£(q, t). We start with the stochastic force £(q, t), the 
noise. It is intuitively evident that such a force always 
exists. In a more mathematical description, this force 
can come either from the environmental influence on the 
network, or from approximations such that the continu- 
ous representation of a discrete process. Here a Gaussian 
and white noise is assumed. In reality, more complicated 
noises, non-Gaussian and colored, can exist. The presen- 
tation formulation in the form of Eq. (4) should provide 
a good starting point. For example, Eq.(4) is already in 
the form in dissipative dynamics™ where colored noises 
have been readily considered. The emphasis on noise in 
the presentation formulation also suggests that metabolic 
fluxes may not be good dynamical variables. The better 
variables may be the numbers of proteins or other macro- 
molecules inside the cell which carry out the metabolic 
task. 

S^q)^. The friction matrix 5(q) and the frictional 
force Sq are self-evident, too. The existence of the fric- 
tion shows that the network has the tendency to ap- 
proach to a steady state. The friction can be the real 
friction in mechanics, or is a representation of the degra- 
dation of proteins or other materials in biology. The 
present dynamical structure theory requires that the fric- 
tion is always associated with the noise according to the 
stochasticity-dissipation relation, Eq.(5). The friction 
and noise are the two opposite sides of stochastic dy- 
namics: the ability to adaptation with friction and the 
ability to optimization with noise. 

T(q)q. The antisymmetric matrix T(q) and the trans- 
verse force Tq are less obvious. In physical sciences, the 
antisymmetric matrix is closely related to quantities in 
physics and chemistry such as the mass, the magnetic 
field, and the rotation. It is a general statement on 



the conjugate relations among state variables. The well- 
known examples of corresponding transverse forces are 
Lorentz and gyroscopic forces. In a complex network, 
the transverse force represents the network ability to re- 
locate the resource from one part of the network to an- 
other. In this sense it is similar to the conversion of en- 
ergy between the kinetic energy and potential energy in 
mechanics, with the aid of a finite mass. The transverse 
force is responsible for oscillatory behaviors in networks. 

u(q). The driving force — du(q) is the gradient of the 
scalar function it(q), the network potential, with respect 
to the network state variable q. The network potential is 
the most important quantity, determining the robustness 
of network dynamics. Knowing this network potential is 
equivalent to knowing the landscape in which the network 
evolves: Its minima determine local equilibrium states, 
and its absolute minimum is the ultimate steady state of 
the network. Therefore important dynamical properties 
of the network can be read from u without explicit solving 
for the time dependent solutions. 

We should point out that the potential w(q) so defined 
is an emerging property of the metabolic network. It has 
no direct thermodynamic meaning as those in physics and 
chemistry. Specifically, it is not the usual free energy 
in thermodynamics. Since the potential it(q) describes 
what the network would eventually like to be under all 
thermodynamic and other constraints, it is a quantity at 
a higher level description than those of free energies to 
describe the network reaction rates. Similar dynamical 
structures as those in Eq.(4) and (5) also exist in other 
branches of biology^. 



V. ILLUSTRATIONS 

In this section we give two illustrations, one mathemat- 
ical and one biological. We first discuss a useful approx- 
imation scheme to find Eq.(4) from Eq.(l). In this way 
the connection between present formulation, Eq. (4) , and 
the classical formation in biochemistry, Eq.(l), becomes 
more clear. Then we brief summarize an successful exam- 
ple of the application the four dynamical element analy- 
sis in an outstanding stability puzzle of a gene regulatory 
network. 



A. First order gradient expansion 

We give an explicit demonstration of how to obtain 
Eq.(4) from Eq.(l) by a so-called gradient expansion to 
the first order derivative of force f with respect to the 
state variables. To be specific, we only consider a two 
dimensional problem. Here q\ and </2 represent numbers 
of two enzymes in a cell. The force in Eq.(l) consists 
of two effects, the production rate and the degradation 
rate: 

/<(q) = /*(q) - /«(q) < = 1, 2 , (10) 
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with the subscripts p and d stand for the production and 
degradation respectively. Under the diffusion approxi- 
mation, the stochastic force isiS*ii 

Ci(q, t) = ,J fi P (<l)ti P (t) + y r J~M)C ld {t) 1 = 1,2, (11) 

with £i p (t), Cid(t) are unity random variables and possible 
correlation among them. Therefore the diffusion matrix 
D can be readily obtained, which is what needed below. 
We remark that the equation similar to above has been 
used in the study of bio-networks 6,9,18 . 

The construction of Eq.(4) from Eq.(l) will be given to 
the lowest order in the gradient expansion. The useful- 
ness of this approximated construction can be illustrated 
for following three reasons. First, in many practical ap- 
plications, this approximation is a good approximation 6 . 
In fact, it is exact in the linear case. Second, the Eq.(7) 
becomes algebraic, instead of a partial differential equa- 
tion. The complete solution of such algebraic equation 
can be founcUS,. Third, several salient features of the 
gauged decomposition becomes apparent without undue 
mathematical complications. An important quantity is 
the force matrix F. According to the definition following 
Eq.(3), 

Fu = dih , F12 = d 2 fi , F 21 = dxh , F 22 = d 2 f 2 . (12) 

Eq.(8) will not change under the gradient approxima- 
tion. In the lowest order gradient approximation, Eq.(7) 
becomes simple. We collect them here: 

/ GF T — FG T = , . 

\G + G T = D ■ { 6) 

In two dimensions the matrix manipulation is particu- 
larly straightforward. The antisymmetric part of the 
auxiliary matrix G = D + Q from Eq.(13) can be found 
to be 

Q = (FD - DF T ) /tr(F) . (14) 
Using the relation 

/ M u M 12 y 1 = 1 f M 22 -M12 \ 
^ M 2 i M 22 J det(M) \ -M 21 M n ) 

The friction matrix S and the antisymmetric matrix T 
can be found according to Eq.(6): 

' u = -I c d4.[G-\q[)f{q[)] 
< S = (^- 2 D ^)/dct(D) . (15) 
^ T = -Q/det(G) 

In two dimensions, det(G) = det(D) + det(Q) = 
D 22 Dn — D\ 2 + Q12 an d is obviously non-negative. 



B. Solving the stability puzzle of a genetic switch 

We have used the scheme here to study the outstanding 
puzzle of stability and efficiency of phage lambda genetic 
switch in a quantitative mannei^. The key was to find 
as accurate as possible the potential function of the gene 
maintenance regulatory network. The dynamical equa- 
tion for such a network can be written down according to 
chemical rate equations. Information at three levels, the 
DNA, protein, and function, needs to be integrated to 
provide a quantitative model. Based on above approach 
we are able to obtain the potential function and achieve a 
quantitative agreement with biological data. We are able 
to make quantitative predictions such as the amount of 
cooperative free energy. We show mathematically that a 
robust and efficient genetic switch can be obtained: po- 
tential function is a measure of robustness and stochastic 
force is responsible for the efficiency. 

We should remark that in reality, both the intrinsic and 
the extrinsic noise coexist. They are equally important 
and measurable^. The stochasticity-dissipation relation 
treats them on the equal footing to determine the gauged 
singular decomposition. We have made used of this fact 
in the modelling of noise in our effort to solve the stability 
puzzle of the genetic switch 6 ^. 



VI. CONNECTIONS TO KNOWN METABOLIC 
NETWORK MODELLING 

In section III we demonstrated that the proposed four 
dynamical component formulation of stochastic differen- 
tial equations is equivalent to the usual chemical rate 
equations. The proposed method offers a direct connec- 
tion to observed data because of the transparent mean- 
ing of the potential and its role in the optimization. In 
the following we briefly discuss its connections to other 
successful modelling methodology in metabolic network 
study. 

1). Flux Balance Analysis (FBA). There are two ingre- 
dients in FBA 21 : The first one is the mass conservation, 
built into the formulation through the stoichiometric con- 
sideration. The second one is a linearizing near a steady 
state state. The first feature is very useful in that it 
eliminates redundant dynamical variables, which would 
be Eq.(l) more tractable. The second feature may be 
viewed as a special case of Eq.(l), by setting the stochas- 
tic force to be zero, and by linearizing the force f around 
its stable fixed point with the given stoichiometric con- 
straints. The justification for the stable fixed point is nat- 
ural because of the existence of heomostasis: The steady 
state of the biological process should be stable under the 
given constraints^. 

Thought it is in general a rather weak constraint, 
FBA has been most easy to implement in practical ap- 
plications, because many matured mathematically tools, 
such as linear programming, can be employed. It is 
currently the working horse in modelling and engineer- 
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ing of metabolic networks, augmented by additional 
considerations 2 ^. 

2) . Metabolic Control Analysis (MCA). To go beyond 
the simple linear case and to study how the network in 
steady state responds to changes in fluxes and its com- 
ponents have been the primary goal of MCA 23 . Proper- 
ties of the architectural structure of the metabolic net- 
work can be revealed by MCA. Undoubtly, MCA may be 
viewed as the study of the structure built into the force f 
in Eq.(l), the situation without dynamics and stochastic 
force. 

3) . Biochemical Systems Theory (BST). A more gen- 
eral and well grounded approach is BST 24 : The force f 
is assumed to polynomial to assist mathematical anal- 
ysis. Time dependent is included into the formulation 
but the stochastic force is typically neglected. Method- 
ologies in cybernetics or control theory are employed in 
BST. This approach is clearly from an engineers point of 
view, providing a great insight for the metabolic network 
engineering. 

4) . Stoichiometric Network Thermodynamic theory 
(SNT). With the concern that there is a tendency to 
ignore the thermodynamic constraints in metabolic net- 
work modelling, the recently developed SNT attempts 
to remedy such shortcomings with explicitly incorpora- 
tion of such constraints 25 . Indeed, because energy and 
entropy play such dominant roles in metabolic network 
dynamics, this is a very desirable progress. Novel results 
have been obtained along this approach^SiSi. Again, the 
starting point of SNT is an equation with the same form 
as of Eq.(l), which shows that our framework is compat- 
ible with SNT, too. It remains to find the connection 
between the potential function in Eq.(4) and those of 
thermodynamic quantities employed in SNT. 

We should point out that in the present article it is im- 
possible to give an adequate survey of vast literature on 
the modelling of metabolic networks. Fortunately, sev- 
eral excellent books already exist 10 ' 21 i 22 i 23 i 24 ' 28 i 29 , where 
ample discussions on FNA, MCA, BST and others can 
be found. The need to incorporate stochastic effects into 



modelling has already been demonstrated by stochastic- 
ity in biological experimental dataSi&S^. 

Another feature should also be pointed here. From 
the present stochastic modelling, Eq.(4) or (1), there are 
obviously two time scales: the very short one character- 
izing the stochastic force £(q, t) or ((q_,t) and the time 
scale on which the smooth functions of potential function 
w(q), degradation (friction) matrix S'(q) and the trans- 
verse matrix T(q) having well defined meaning. This cor- 
responds nicely to the hierarchical structure abundant in 
metabolic pathway analysis 31 . 



VII. CONCLUSION 

In this paper we have proposed a general systems ap- 
proach to model metabolic network dynamics. It is based 
on a special form of stochastic differential equations. We 
have demonstrated that it is equivalent to the classical 
chemical rate equations approach with noise. The con- 
nections to various existing modelling methodologies are 
pointed out: In each of their valid description regions 
our method is compatible. Hence it may provide a uni- 
fying mathematical framework with the particular use- 
ful potential function. Its usefulness has already been 
demonstrated in the study of a gene regulatory network. 
It remains to be demonstrated that additional biologi- 
cal insights in metabolic network study can be obtained. 
This is the task been undertaking. 
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